Traffic of interacting ribosomes on mRNA during protein synthesis: 
effects of chemo-mechanics of individual ribosomes 
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Many ribosomes simultaneously move on the same messenger RNA (mRNA), each synthesizing 
separately a copy of the same protein. In contrast to the earlier models, here we develop a "unified" 
theoretical model that not only incorporates the mutual exclusions of the interacting ribosomes, but 
also describes explicitly the mechano-chemistry of each of these macromolecular machines during 
protein synthesis. Using analytical and numerical techniques of non-equilibrium statistical mechan- 
ics, we analyze the rates of protein synthesis and the spatio-temporal organization of the ribosomes 
in this model. We also predict how these properties would change with the changes in the rates 
of the various chemo-mechanical processes in each ribosome. Finally, we illustrate the power of 
this model by making experimentally testable predictions on the rates of protein synthesis and the 
density profiles of the ribosomes on some mRNAs in E-coli. 
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I. INTRODUCTION 



Synthesis of each protein from the corresponding mes- 
senger RNA (mRNA) is carried out by a ribosome [1] 
and the process is referred to as translation (of genetic 
code). Ribosome is one of the largest and most sophisti- 
cated macromolecular machines within the cell J2, l3| . It 
is not merely a "protein-making motor protein" [J, |5| but 
it also serves as a "workshop" which provides a platform 
where a coordinated action of many tools take place for 
the synthesis of each of the proteins. What makes it very 
interesting from the perspective of statistical physics is 
that most often many ribosomes move simultaneously on 
a single mRNA strand while each synthesizes a separate 
copy of the same protein. Such a collective movement 
of the ribosomes on a single mRNA strand has superfi- 
cial similarities with vehicular traffic [gland, therefore, 
will be referred to as ribosome traff ic 3 , [ a , M ■ All the 
earlier works [H El El G3 EI EE EllHEi treat ri- 
bosome traffic as a problem of non-equilibrium statistical 
mechanics of a system of interacting "self-driven" hard 
rods. But, strictly speaking, a ribosome is neither a par- 
ticle nor a hard rod; its mechanical movement along the 
mRNA track is coupled to its internal mcchano-chemical 
processes which drive the synthesis of the protein. Thus, 
one serious limitation of the earlier models is that these 
cannot account for the effects of the intra-ribosome chem- 
ical and conformational transitions on their collective 
spatio-temporal organization. 

In this paper we develop a model that not only in- 
corporates the inter-ribosome steric interactions (mutual 
exclusion), but also captures explicitly the essential steps 
in the intra-ribosome chemo-mechanical processes. From 
a physicist's perspective, our model is a biologically mo- 
tivated extension of some earlier models developed for 
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studying the collective spatio-temporal organization in a 
non-equilibrium system of interacting "self-driven" hard 
rods. However, in contrast to the earlier models, each of 
the rods in our model has several "internal" states which 
capture the different chemical and conformational states 
of an individual ribosome during its biochemical cycle. 

Our modelling strategy for incorporating the biochemi- 
cal cycle of the ribosomes is similar to that followed in the 
recent work [l9j] on single-headed kinesin motors KIF1A. 
However, the implementation of the strategy is more dif- 
ficult here not only because of the higher complexity 
of composition, architecture and mechano-chemical pro- 
cesses of the ribosomal machinery but also because of the 
sequence heterogeneity of the mRNA track [2^, [21] . Our 
approach is based on a stochastic chemical kinetic model 
that describes the dynamics in terms of master equa- 
tions. But our model makes no commitments to either 
power stroke or Brownian ratchet mechanism [IE HE H3| 
of molecular motors. 

The paper is organized as follows: Because of the in- 
terdisciplinary nature of the topic investigated in this 
paper, we present in section [IT] a summary of the essen- 
tial biochemical and mechanical processes during a com- 
plete operational cycle of a single ribosome. We present 
a brief critical review of the earlier models in section IIIII 
followed by a description of our model in section [IV] so as 
to highlight the novel features of our model. We report 
our results on this model with periodic boundary condi- 
tions (PBC) in section [V] and those with open boundary 
conditions (OBC) in section IVT1 Wc summarize the main 
results and draw conclusions in section IVTT1 



II. SUMMARY OF THE ESSENTIAL 
CHEMO-MECHANICAL PROCESSES 

A protein is a linear polymer of amino acids which 
are linked together by peptide bonds and, therefore, of- 
ten referred to as a polypeptide. An mRNA is a linear 
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FIG. 1: A pictoral depiction of three major steps in the 
chemo-mechanical cycle of a single ribosome. The larger and 
smaller subunits have been depicted as two rectangles. The 
complementary shapes of the vertical tips and dips merely 
emphasize the codon-anticodon matching. 



polymer of nucleotides and triplets of nucleotides form 
one single codon. The stretch of mRNA between a start 
codon and a stop codon serves as a template for the syn- 
thesis of a polypeptide. The process of translation itself 
can be divided into three main stages: (a) initiation, dur- 
ing which the ribosomal subunits assemble on the start 
codon on the mRNA strand, (b) elongation, during which 
the nascent polypeptide gets elongated by the formation 
of peptide bonds with new amino acids, and (c) termina- 
tion, during which the process of translation gets termi- 
nated at the stop codon and the polypeptide is released. 
Initiation or termination can be the rate-limiting stage 
in the synthesis of a protein from the mRNA template 
[lij . In this paper we shall be concerned mostly with the 
process of elongation. 

The specific sequence of amino acids on a polypeptide 
is dictated by the sequence of codons on the correspond- 
ing template mRNA. The dictionary of translation relates 
each type of possible codons with one of the 20 species of 
amino acids; these correspondences are implemented by 
a class of adaptor molecules called tRNA. One end of a 
tRNA molecule consists of an anti-codon (a triplet of nu- 
cleotides) while the other end carries the cognate amino 
acid (i.e., the amino acid that corresponds to its anti- 
codon). Since each species of anti-codon is exactly com- 
plimentary to a particular species of codon, each codon 
on the mRNA gets translated into a particular species of 
amino acid on the polypeptide. A tRNA molecule bound 
to its cognate amino acid is called aminoacyl-tRNA (aa- 
tRNA). 
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FIG. 2: A schematic representation of the biochemical cycle 
of a single ribosome during the elongation stage of transla- 
tion in our model. Each box represents a distinct state of the 
ribosome. The index below the box labels the codon on the 
mRNA with which the smaller subunit of the ribosome binds. 
The number above the box labels the biochemical state of 
the ribosome. Within each box, 1(0) represents presence (ab- 
sence) of tRNA on binding sites E, P, A respectively. 1* is 
a elongation factor (EF)-Tu bound tRNA and G is a EF-G 
GTPase. The symbols accompanied by the arrows define the 
rate constants for the transitions from one biochemical state 
to another. As explained in section ITVl the dashed arrow rep- 
resents the approximate pathway we have considered in the 
simplified dynamics of our model. 



Each ribosome consists of two parts which are usually 
referred to as the larger and the smaller subunits. There 
are four binding sites on each ribosome. Of these, three 
sites (called E, P, A), which are located in the larger 
subunit, bind to aminoacyl-tRNA (aa-tRNA), while the 
fourth binding site, which is located on the smaller sub- 
unit, binds to the mRNA strand. The translocation of 
the smaller subunit of each ribosome on the mRNA track 
is coupled to the biochemical processes occuring in the 
larger subunit. 

Three major steps in the biochemical cycle of a ribo- 
some are sketched schematically in figHJ In the first, the 
ribosome selects an aa-tRNA whose anticodon is exactly 
complementary to the codon on the mRNA. Next, it cat- 
alyzes the reaction responsible for the formation of the 
peptide bond between the existing polypeptide and the 
newly recruited amino acid resulting in the elongation 
of the polypeptide. Finally, it completes the mechano- 
chemical cycle by translocating itself completely to the 
next codon and is ready to begin the next cycle. Elon- 
gation factors (EF), which are themselves proteins, play 
important roles in the control of these major steps (see 
the fig[l|) which require proper communication and coor- 
dination between the two subunits. Moreover, some spe- 
cific steps in the mechano-chemical cycle of a ribosome 
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are driven by the hydrolysis of guanosine triphosphate 
(GTP) to guanosine diphosphate (GDP). 

The detailed chemo-mechanical cycle of a ribosome is 
drawn schematically in figH] Let us begin the biochem- 
ical cycle of a ribosome in the elongation phase with 
state 1, (figure \2§ where a tRNA is bound to the site 
P. A tRNA-EF-Tu complex (a macromolecular complex 
of a tRNA and an elongation factor Tu) now binds to 
site A and, after the correct recognition of the cognate 
aa-tRNA through proper codon-anticodon matching, the 
system makes a transition from the state 1 to the state 
2. As long as the EF-Tu is attached to the tRNA, codon- 
anticodon binding can take place, but the peptide bond 
formation is not possible. The EF-Tu has a GTP part 
which is hydrolyzed to GDP, driving the transition from 
state 2 to 3. Following this, a phosphate group leaves, 
resulting in the intermediate state 4. This hydrolysis, fi- 
nally, releases the EF-Tu, and then the peptide bond for- 
mation becomes possible. After this step, another elon- 
gation factor, namely, EF-G, in the GTP bound form, 
comes in contact with the ribosome. This causes the tR- 
NAs to shift from site P to E and from site A to P, site 
A being occupied by the EF-G, resulting in the state 5. 
Hydrolysis of the GTP to GDP then releases the EF-G 
and this is accompanied by the transition of the system 
from the state 5 to the state 6. The transition from the 
state 6 to the state 7 is accompanied by conformational 
changes which are responsible for the forward movement 
of the smaller subunit by one step. Finally, the tRNA 
on site A is released, resulting in the completion of one 
biochemical cycle; in the process the ribosome completes 
its forward movement by one codon (i.e., one step on the 
lattice). In each cycle of the ribosome during elongation, 
the search for the cognate tRNA is the rate limiting step 
[25l [2^|. the corresponding rate constant being cu a . 



III. BRIEF REVIEW OF THE EARLIER 
MODELS 

To our knowledge, MacDonald, Gibbs and coworkers 
[lrllTiT] developed the first quantitative theory of simulta- 
neous protein synthesis by many ribosomes on the same 
mRNA strand. The sequence of codons on a given mRNA 
was represented by the corresponding sequence of the eq- 
uispaced sites of a regular one-dimensional lattice. The 
details of molecular composition and architecture of the 
ribosomes was ignored in this model. Instead, each of the 
ribosomes was modelled by a "self-propelled particle" of 
size £ in units of the lattice constant; thus, £ is an integer. 
On the lattice the steric interaction of the ribosomes was 
taken into account by imposing the condition of mutual 
exclusion, i.e., no site of the lattice can be simultaneously 
covered by more than one particle. 

The dynamics of the system was formulated in terms 
of the following update rules: 

An extended particle (effectively, a hard rod), whose for- 
ward edge is located at the site i, can hop forward by 



one lattice spacing with the forward hopping rate 
provided the target site is not already covered another 
extended particle. Moreover, initiation and termination 
were assumed to take place with the corresponding rates 
a and /3, respectively, which are not necessarily equal to 
any of the other rate constants. 

For the sake of simplicity of analytical calculations, 
one usually replaces this intrinsically inhomogeneous pro- 
cess by a hypothetical homogeneous one by assuming 
[ToL [ll| that q l = q, irrespective of i. In such spe- 
cial situations, this model reduced to the totally asym- 
metric simple exclusion process (TASEP) without any 
defect or disorder (27[, except that the allowed size of 
the "extended" hopping particles is a multiple of the 
lattice spacing. TASEP is one of the simplest mod- 
els of systems of interacting driven particles [28[. Note 
that these TASEP-like earlier models of ribosome traffic 
[lEllllMIIilMIMIMSli! capture the effects of 
all the chemical reactions and conformational changes, 
which lead to the translocation of a ribosome from one 
codon on the mRNA to the next, by the single parameter 

The steady-state flux J of the ribosomes is defined as 
the average number of the ribosomes crossing a specific 
codon (selected arbitrarily) per unit time. Because of 
the close analogy with vehicular traffic, we shall refer to 
the flux-density relation as the fundamental diagram [|| . 
In the context of ribosomal traffic, the position, average 
speed and flux of ribosomes have interesting interpre- 
tations in terms of protein synthesis. The position of 
a ribosome on the mRNA also gives the length of the 
nascent polypeptide it has already synthesized. The av- 
erage speed of a ribosome is also a measure of the av- 
erage rate of polypeptide elongation. The flux of the 
ribosomes gives the total rate of polypeptide synthesis 
from the mRNA strand, i.e., the number of completely 
synthesized polypeptides per unit time. 

The rate of protein synthesis and the ribosome density 
profile in the model developed by Macdonald et al. [ToL [Tl1| 
as well as in some other closely related models have been 
investigated in detail. PBC are less realistic than OBC 
for capturing protein synthesis by a theoretical model. 
Nevertheless, if one imposes PBC on this simplified ver- 
sion of the model, the steady-state flux of the ribosomes 
is given by [l(| Ha] 



J = q 



P(l -Pi) 
l-p(£-l) 



(1) 



where p is the number density of the ribosomes; if N is 
the total number of ribosomes on the lattice of length 
L, then p = N/L. The corresponding average speed of 
the ribosomes is given by < v >— J/ p. In the special 
case £ = 1 the expression (JTJ) reduces to the well known 
formula 



J = Q Pi 1 - P) 



(2) 



for the steady-state flux in the TASEP. Comparing equa- 
tion (JTJ) with @, p/[l —p(£— 1)] has often been identified 
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(37j as an effective particle density while the correspond- 
ing effective hole density is given by 1 — p£. The cor- 
responding phenomenological hydrodynamic theory [l3[ 
has also been derived [H, [39[ from the TASEP-like dy- 
namics of the hard rods of size £ on the discrete lattice. 

The fundamental diagram implied by the expression 
([T]) exhibits a maximum at the density p m and the value 
of flux at this maximum is J m where 

Pm = — F 7= an( i Jm = r=~ • (3) 

Only in the special case £ = 1, this fundamental dia- 
gram is symmetric about p = 1/2; the maximum shifts 
to higher density with increasing I. 



IV. THE MODEL 



(a) 




B 

FIG. 3: A schematic representation of the model, (a) A car- 
toon of a single ribosome that explicitly shows the three bind- 
ing sites E, P and A on the larger subunit which is represented 
by the ellipsoidal lobe. The rectangular lower part repre- 
sents the smaller subunit of the ribosome. (b) The mRNA is 
represented as a one-dimensional lattice where each site cor- 
responds to one single codon. The smaller subunit of each 
ribosome covers I codons (I = 2 in this figure) at a time. 

Our model is shown schematically in fig[3] We rep- 
resent the single-stranded mRNA chain, consisting of L 
codons, by a one-dimensional lattice of length L + £ — 1 
where each of the first L sites from the left represents 
a single codon. We label the sites of the lattice by the 
integer index i; the site i = 1 represents the start codon 
while the site i — L corresponds to the stop codon. 
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25 


25 


0.0028 


10 


10 


2.4 



TABLE I: Rate constants obtained from experimental data 
for E-coli 53, EH- 



Our model differs from all earlier models in the way 
we capture the structure, biochemical cycle and translo- 
cation of each ribosome. The small sub-unit of the ri- 
bosome, which is known to bind to the mRNA, is rep- 
resented by an extended particle of length £ which is ex- 
pressed in the units of the size of a codon (see figE])- 
Thus, in our model, the small subunit of each ribosome 
covers £ codons at a time; no lattice site is allowed to 
be covered simultaneously by more than one overlapping 
ribosome. Irrespective of the length £, each ribosome 
moves forward by only one site in each step as it must 
translate successive codons one by one. 

Since our model is not intended to describe initiation 
and termination in detail, we represent these two pro- 
cesses by the parameters a and (3 respectively. Whenever 
the first £ sites on the mRNA are vacant this group of 
sites is allowed to be covered by a fresh ribosome with 
the probability a in the time interval At (in all our nu- 
merical calculations we take At — 0.001 s). Similarly, 
whenever the rightmost £ sites of the mRNA lattice are 
covered by a ribosome, i.e., the ribosome is bound to the 
L-th codon, the ribosome gets detached from the mRNA 
with probability (3 in the time interval At. Since a is 
the probability of attachment in time At, the probability 
of attachment per unit time (which we call co a ) is the 
solution of the equation a = 1 — e -"° xAt ( se e appendix 
A for the detailed explanation) . Similarly, we also define 
u)p which is the probability of detachment per unit time. 

In the elongation stage, we have identified seven major 
distinct states of the ribosome in each cycle which have 
been described in detail in section [II] (shown schemati- 
cally in figUl)- However, in setting up the equations be- 
low, we further simplify the model. Throughout this pa- 
per, we replace the pathway 5^6— » 7 — > 1 by an effec- 
tively direct single transition 5 — > 1, with rate constant 
ujh2 (shown by the dashed line in figHJ. This simplifica- 
tion is justified by the fact that the transitions 5 — > 6 and 
6^7 are purely "internal" , and do not seem to depend 
on the availability of external molecules like elongation 
factors, or OTP or aa-tRNA. 

The typical values of the rate constants have been 
extracted from empirical data for the bacteria E-coli 
(40l l4ll ]. Moreover, since there is no significant differ- 
ence in the structures of the two elongation factors and 
since their binding mechanisms are also similar Q, we 
assume that the rate constants Uhi and coh2 are equal. 
The values of the rate constants used in our calculations 
are listed in table HI 

The lifetime of a typical eukaryotic mRNA is of the or- 
der of hours whereas the time taken to synthesize an en- 
tire protein by translating the mRNA is of the order of a 
few minutes. Consequently, most often protein synthesis 
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takes place under steady-sate conditions. Therefore, al- 
though we shall formulate time-dependent equations for 
protein synthesis, we shall almost exclusively focus on 
the steady-state properties of these models in this paper. 

Most of our analytical calculations have been per- 
formed in the mean-field approximation. In order to test 
the accuracy of the approximate analytical results, we 
have also carried out computer simulations of our model. 
Since we found very little difference in the results for 
systems size L = 300 and those for larger systems, all 
of our production runs were carried out using L = 300. 
We have used random sequential updating which corre- 
sponds to the master equations formulated for the analyt- 
ical description. In each run of the computer simulations 
the data for the first five million time steps were dis- 
carded to ensure that the system, indeed, reached steady 
state. The data were collected in the steady state over 
the next five millon time steps. Thus, each simulation 
run extended over a total of ten million time steps. For 
example, the average steady-state flux was obtained by 
averaging over the last five million time steps. 



V. RESULTS UNDER PERIODIC BOUNDARY 
CONDITIONS 

Typically, a single ribosome itself covers about twelve 
codons (i.e., I = 12), and interacts with others by mutual 
exclusion. The position of such a ribosome will be referred 
to by the integer index of the lattice site covered by the 
leftmost site of the smaller subunit. Main results for the 
special case £ = 1 are given separately in the appendix 
B. 



A. Mean field analysis under periodic boundary 
conditions 

Let Pn(i) be the probability of finding a ribosome at 
site i, in the chemical state fi. Also, P(i) — J2u,=i 
is the probability of finding a ribosome at site i, in any 
state. Let P(i\j) be the conditional probability that, 
given a ribosome at site i, there is another ribosome at 
site j. Then, Q(i\j) = 1 — P(i\j) is the conditional prob- 
ability that, given a ribosome in site i, site j is empty. In 
the mean-field approximation, the Master equations for 
the probabilities P^(i) are given by 



dPi(i) 
dt 



UJh2P 5 (i - l)Q(kl|i - 1 +£ ) + U p P 2 (i) - LO a Pl (i) 

(4) 



dP 2 (i) 
dt 



- v a Pi(i) ~ ( w p +w/ l i)P 2 (i) (5) 



dt 



u hl P 2 {i) - k 2 P 3 (i) 



(6) 



dP A {i) 
dt 



= k 2 P 3 {i) - u g P 4 (i) 



dP 5 (i) 
dt 



oj g P 4 (i) - Lu h2 P 5 (i)Q(i\i + £) 



(7) 



(8) 



Note that not all of the five equations (I!])-© are inde- 
pendent of each other because of the condition 



P(i) 



N 

T 



(9) 



In our calculations below, we have used the equations 
©-([H]) as the five independent equations. Mean field ap- 
proximation has entered through our implicit assumption 
that the probability of there being a ribosome at site i 
is not affected by the presence or absence of other ribo- 
somes at other sites. 



B. Steady state properties under periodic 
boundary conditions 

In the steady state, all the become independent 

of time. Moreover, if PBC are imposed (i.e., the lattice 
effectively forms a ring), no site has any special status 
and the index i can be dropped. The corresponding flux 
of the ribosomes J can then be obtained from 



J = Lo h2 P 5 Q(i\i + £). 



(10) 



using the steady-state expressions for Q(i\i + £) and P5. 

Because of the translational invariance in the steady 
state, we have Q(i\j) = Q(l\j — i + 1). Therefore, we 
now calculate Q(l|l + I): given a ribosome at the site 
2 = 1, what is the probability that the site i = £ + 1 is 
empty? Since it is given that the site i = 1 is occupied 
by a ribosome, the remaining N — 1 ribosomes must be 
distributed among the remaining L — £ sites. Let us in- 
troduce the symbol Z(L, N, £) to denote the number of 
ways of arranging the N ribosomes and L — N£ gaps. 
Obviously, 



Z{L,N,i) 



(N + L-Nl)\ 
N\{L-Nl)\ 



(11) 



In case it is given that one ribosome occupies i = 1, the 
total number of configurations would be Z(L—£, N—l, £). 
Of these, we wish to find the number of those configu- 
rations where i = £ is also occupied; this is given by 
Z{L — 2£,N— 2, £). Therefore, the probability that i = £ 
is occupied, given that i = 1 is occupied, is given by 



P(l|*+1) = 



Z(L - 21, N - 2, 1) 
Z{L- i,N -1,1) 

N - 1 
L + N - N£ - 1 



(12) 
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Hence, 



Q(i\i + £) 



L-Nl 
L + N — N£ — 1 ' 



(13) 



Solving the equations (HE]) in the steady state under 
PBC, we get 



Pa = 



1 + 



Q. h2 (L-Nl) 
L+N-Nl-l 



where, 



with 



*eff 



hi 



1 1 

Ulg k 2 



Uh2/k e ff- 



1 



(14) 



(15) 



(16) 



Note that is an effective time that incorporates the 
delays induced by the intermediate biochemical steps in 
between two successive hoppings of the ribosome from 
one codon to the next. Therefore, k e ff — > oo implies 
short-circuiting the entire biochemical pathway so that a 
newly arrived ribosome at a given site is instantaneously 
ready for hopping onto the next site with the effective 
rate constant LO) l2 . 

Using expressions (fT3|) and (fl~4|) in (fTUf and the defini- 
tion p = N I L for the number density, we get 



J = 



(1 



+ n h2 (i- P z) 



(17) 



Note that J vanishes at p = and for all p > p max = 1/1 
because at the density p max the entire mRNA in fully 
covered by ribosomes. Moreover, in the special situation 
where k e ff — > oo, but ojh 2 = Q remains non-zero and 
finite, Slh2 — > and the expression (fT7|) reduces to the 
expression ([T]). 

The flux obtained from (fT7|) has been plotted in figure 
(|4]) for two values of I. This trend of variation with £ was 
also observed in the pioneering work of MacDonald et al. 
[lol ] in their simpler model. By differentiating equation 
([TT)) . we obtain that the density p* corresponding to the 
maximum of the flux is the solution of the equation: 

p 2 £(l - £ - Q h2 £) + 2p£(l + Q h2 ) - (1 + fl h2 ) - (18) 



and, hence, 




y/£(i + n h2 ) + 1 



(19) 



We recover equation (|3|) for p m in the appropriate limit 
^/i2 — * 0. Our theoretical predictions in fig[5]are also in 
good agreement with the corresponding simulation data. 

In order to see the effects of varying the rates of some 
of the biochemical transitions, we have plotted the fun- 
damental diagram in figH] for two different situations, 
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p (number per codon) 



FIG. 4: Flux of ribosomes with I — 3, 12, under periodic 
boundary conditions, plotted against density. The curves cor- 
respond to the analytical expressions obtained in the mean- 
field (MF) approximation whereas the discrete data points 
have been obtained by carrying out computer simulations 
(Sim). Values of all the parameters, except £, are same as 
those listed in table I. 
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FIG. 5: Flux of ribosomes with i — 12, under periodic bound- 
ary conditions, plotted against density. The curves corre- 
spond to the analytical expressions obtained in the mean-field 
(MF) approximation whereas the discrete data points have 
been obtained by carrying out computer simulations (Sim). 
Values of all the parameters, which are not shown explicitly 
on the figure, are identical to those in table I. 



namely, u>hi = 10uJh 2 with ujhi = 10 s _1 and Wh 2 = lOa^i 
with u>hi = 10 s _1 . The fundamental diagrams in these 
two situations turn out to be almost identical; this is 
a consequence of the fact that for the set of parameter 
rages used in this figure, neither u>hi nor ujh 2 corresponds 
to the rate limiting process. 

We have plotted the fundamental diagrams of the 
model in fig [S] for three different different values of u) a , 
namely, u a — 2.5 s^ 1 , uo a — 25 s _1 and uj a = 250 
s _1 using both mean-field theory and computer simu- 
lations. The results show that at sufficiently small values 
of u> a , where the availability of the tRNA is the rate- 
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FIG. 6: Same as in fig[5] except that different curves corre- 
spond to different values of w a . Values of all the parameters, 
which are not shown explicitly on the figure, are identical to 
those in table I. 



FIG. 8: Same as in fig[5] except that different curves corre- 
spond to different values of &2- Values of all the parameters, 
which are not shown explicitly on the figure, are identical to 
those in table I. 



VI. RESULTS UNDER OPEN BOUNDARY 
CONDITIONS 
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FIG. 7: Same as in fig[5] except that different curves corre- 
spond to different values of ui g . Values of all the parameters, 
which are not shown explicitly on the figure, are identical to 
those in table I. 



limiting process, the flux increases rapidly with increas- 
ing ui a . However, the rate of this increase decreases with 
increasing uj a and, eventually, the flux essentially satu- 
rates. This saturation of flux occurs when Lu a is so large 
that the availability of tRNA is no longer the rate limit- 
ing process. A similar trend of variation of flux with w g 
is observed in figlT] when tu g is varied over three orders 
of magnitude. In contrast, the flux has been observed to 
vary at a significant rate even at the highest values of 
&2, when it is varied over three orders of magnitude (see 
figO indicating that saturation of flux with respect to k 2 
variation sets in at even higher values of k 2 ). 



An OBC is more realistic than a PBC for describing ri- 
bosome traffic on mRNA. The parameters a and /3, which 
are associated, respectively, with initiation and termina- 
tion of translation play significant roles in the system 
under OBC. 



Mean field analysis under open boundary 
conditions 



In this subsection we calculate the flux of ribosomes 
(and, hence, the rate of protein synthesis) using a mean 
field theoretical approach similar to that developed by 
Shaw et al. [l4| . The approximation involved is that the 
conditional probability of site i + £ being empty, given 
that site i has a ribosome in it, is replaced simply by the 
probability of site i being empty, given no other condi- 
tion. If P(i) is the probability of there being a ribosome 
at site i, then the probability of there being a hole at site 

j is given by 1 — X)s=o P(J ~ s )- ^ ^ s now straightforward 
to set up the master equations for the probabilities P^,(i) 
in the mean-field approximation: 



dPi(l) 
dt 



= u a (l - V P(s)) + u p P 2 (l) - w„Pi(l) (20) 



dPi(i) = m h 2P 5 (i - 1)(1 - ELi p ( l ~ 1 + s )) 
dt 1 - £^ =1 P(i - 1 + s) + P{i -1+t) 

+U! p P 2 {l) - u a P x {i) (21) 

(Mi) 



dP 2 (i) 
dt 



u a P 1 (i)-(Lo p +Lo hl )P 2 (i) (22) 
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dt 
dt 



= w M P 2 (i)-fc 2 P 3 (i) 

= k 2 P 3 (i) - Ul g P 4 (l) 



(23) 
(24) 



0.12 



dP 5 (i) 
dt 



UgPiW - f — (25) 

l-Y,LiP{i + s) + P{i + l) 

w g P A (N) - ujpP^N) (26) 



dP 5 (N) 
dt 



B. Steady state properties under open boundary 
conditions 

The flux is given by J — uj a {\ - £f =n P(s)). This 
flux has been computed numerically by solving equations 
(|20j|26[) ; the results are shown by the continuous curves 
in figures [Ha) and[9]^b). These mean-field estimates are 
in excellent agreement with the corresponding numerical 
data obtained from computer simulations of the model. 
Moreover, the rates of protein synthesis corresponding to 
the typical rate constants given in table I are in the same 
order of magnitude as those observed experimentally [l[ . 
The average density profiles observed at several values of 
uj a and ojh are also shown in the insets of figs[9]Ja) and 
(b), respectively. 

Figures [9^a) and (b) show how the current gradually 
increases and, finally, saturates as u> a (in (a)) and u>h (in 
(b)) increases; the saturation value of the current is nu- 
merically equal to the maximum current obtained in the 
corresponding closed system with periodic boundary con- 
ditions. The average density profiles in the insets of the 
figures [9ja) and (b) establish that the average density 
of the ribosomes increases with increasing u> a , but de- 
creases with increasing u>h, gradually saturating in both 
the cases. These observations are consistent with the 
scenario of phase transition from one dynamical phase to 
another, as predicted by the extremal current hypothesis 
which will be considered later in this section. 



C. Effects of inhomogeneity of the mRNA track of 
real gene sequences 

In a real mRNA the nucleotide sequence is, in gen- 
eral, inhomogeneous. First of all, different codons appear 
on an mRNA with different frequencies. Moreover, in a 
given cell, not all the tRNA species, which correspond to 
different codon species, are equally abundant. Interest- 
ingly, because of evolutionary adaptations, the concen- 
trations of tRNA species which correspond to rare codons 
are also proportionately low [42]. Thus, sequence inho- 
mogeneity on a real mRNA can have important effects 
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FIG. 9: Flux of ribosomes plotted against (a) uj a and (b) 
ujh for the genes err (170 codons) and cysK (324 codons) of 
Escherichia coli K-12 strain MG1655, as well as the cor- 
responding curve for a homogenous mRNA strand of 300 
codons. The insets show the average density profiles on a 
hypothetical homogeneous mRNA track for four different val- 
ues of (a) uj a and (b) U)h, for fixed tu a = 25 s _1 . 



We now extend our homogeneous model assuming that 
the rate constant oj a of the attachment of the tRNA to 
the site A of the ribosome is site-dependent (i.e., depen- 
dent on the codon species). More precisely, for a ribo- 
some located at the i-th site, we multiply the numerical 
value of u) a , which we used earlier for the hypothetical 
homogeneous mRNA, by a multiplicative factor that is 
proportional to the relative concentration of the tRNA 
associated with the z-th codon [42l |44]|. 

A lot of work on TASEP with quenched random space- 
dependent hopping rates [H, [H, [4^, [H, [H, H3, j5j, [53, f5^| 
and Brownian motors with quenched disorder [5J, l55l l56l . 
[o7j has been reported. Similarly, effects of randomness 
on some stochastic chemical kinetic models of molecu- 
lar motors have also been investigated [13, Hl|. However, 
the nucleotide sequence on any real DNA or mRNA is not 
random. But, to our knowledge, for the inhomogeneous, 
but correlated, gene sequences no analytical technique is 
available at present for the calculation of the speed of the 
associated molecular motors. For example, all the theo- 
retical schemes developed so far for RNA polymerase mo- 
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a = p 




P=l-P + 



the extrema principle states that 




P + 



FIG. 10: Incorporating a and /3 through two reservoirs with 
appropriate densities. 



tors [58|, [59( by taking into account the actual sequence of 
the specific DNA track, have to be implemented numer- 
ically. Even in the context of earlier TASEP-like mod- 
els of protein synthesis, almost all the theoretical results 
on the effects of sequence inhomogeneities have been ob- 
tained by computer simulations [1 SI ] - Therefore, we have 
been able to study the effects of sequence inhomogeneities 
of real codon sequences on the rate of protein synthesis 
in our model only numerically by carrying out computer 
simulations. 

In our numerical studies, we focus on genes of Es- 
cherichia coli K-12 strain MG1655 [(5(J. We use the hy- 
pothesis mentioned above for the choice of the numerical 
values for the different species of codons. The results of 
our computer simulations are plotted with discrete points 
in figlU The lower flux observed for real genes, as com- 
pared to that for a homogeneous mRNA, is caused by 
the codon specificity of the available tRNA molecules. 



D. Phase diagrams from extremal current 
hypothesis 



We shall treat a, u> a , u>hi and u>h2 as the experimen- 
tally controllable parameters. We shall plot the phase di- 
agrams of the model in planes spanned by pairs of these 
parameters. We shall plot these phase diagrams using 
equation (jTTJ) , and an extremum principle which was orig- 
inally introduced by Krug [6l| , stated in its general form 
by Popkov and Schiitz [62| and effectively utilized in sev- 
eral later works [63], [§4| in the context of driven diffusive 
lattice gas models. 

In this approach, one imagines that the left and right 
ends of the system are connected to two reservoirs with 
the appropriate number densities p_ and p + of particles 
(ribosomes) so that, assuming the same jumping rates as 
in the bulk, the rates a and j3 are incorporated into the 
model (see fig HO]). 

The extrema principle then relates the flux J in the 
open system to the flux J{p) for the corresponding closed 
system (i.e., the system with periodic boundary condi- 
tions) with the same dynamics. In the limit L — > 00 [62l |. 



J 



max J(p) if p- > p > p+ 
rain J(p) if p- < p < p+ 



(27) 



In the present context of our model the expression ([171) 
for J(p) exhibits a single maximum at p = p* where p* 
is given by the equation (|19| . Moreover, we take p+ = 0, 
i.e., (3 = 1, because we assume that the ribosome is re- 
leased from the mRNA as soon as it reaches the stop 
codon; this is justified by the fact that, normally, ter- 
mination is not the rate limiting step in the process of 
protein synthesis. Therefore, the extremal current hy- 
pothesis implies that in our model 



J = max J(p) if P- > p* 



(28) 



All the results derived in this section exploiting this ex- 
tremum principle are approximate because the expression 
(jTTJ) , which we use for the expression of J(p), has been 
derived in the mean-field approximation. Next, following 
arguments similar to those followed in all earlier appli- 
cations of this extremum principle, we now derive the 
appropriate expressions for p_ . 

Consider a closed system with L sites. Given that a 
sequence of £ successive sites are empty, the total number 
of ways in which N ribosomes can be distributed over 
the remaining L — £ sites is simply Z(L — £,N,£). Of 
these, Z(L — 2£, N — 1, £) configuations have a ribosome 
in the adjacent £ sites to the left of these empty £ sites. 
Let us use the symbol P(11...1| 0Q^Q ) for the conditional 

probability that, given a sequence of £ successive empty 
sites, there will be a ribosome in the adjacent £ sites to 
its left. Then, from the above considerations, 



Pdl-.-ll OO.-O ) 



Z(L-2£,N -1, 
Z{L-£,N, £) 



N 



L — £ + N - N£ 



In terms of the number density p, this probability can 
expressed as 



Pfll-ll OO-O ) = 



p+l-f-pt 



Moreover, solving the equations (|20II26[) we find 

1 



1 + fi 



1:2 



(30) 



(31) 



Therefore, if pi um v i s the probability that, given a se- 
quence of £ successive empty sites, a ribosome will hop 
onto it in the next time step, we have 

pjump = pf 11 ,.. 1 | oo„.Q ) x P 5 x uj h2 (At) (32) 



where P(1L ; _1| 00...0 ) and P 5 are given by ([3U| and ^T\ , 
respectively. 
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FIG. 11: Phase diagram in a — P u}a plane. Pu Ja is the proba- 
bility of attachment of a tRNA in time At — 0.001 s, and is 
related to uj a by equation ((36} . This diagram has been plotted 
for uj h \ = Lu h2 = uj h = 10 sec -1 . 
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FIG. 12: Phase diagram in P UJh — P^ a plane. P^ h is the 
probability of hydrolysis in time At = 0.001 s, and is related 
to ujh by equation (136[) . 



Now, going back to the open system, p_ is the solution 
of the equation a = pi um P and, hence, we get 



9- = 



a(i-£)(i + n M ) 



Pu 



a(l+n M )(l-i) 



(33) 



where P Uh is the probability of hydrolysis in time At. 
In the special case fc e // — > oo while LOfa — q remains 
finite and non-zero, i.e., flh2 — * 0, p- 



l+«-l)a » Which 

is identical to the expression derived earlier by Lakatos 
and Chou [l2j for this special case. 

Thus, the boundaries between various phases have 
been obtained by solving the equation 

p~(a, u a , LOhUUhi) = 9*(a,u a ,Uhi,u)h2) (34) 

numerically using p* and p_ obtained, respectively, from 
equations ([IT?]) and (|33[) ; two typical phase diagrams have 
been plotted in in figures [11] and assuming [J, 0, Hl| 



dj/ji = u>h2 = ^>h- Each of these phase diagrams show 
two phases namely, a maximal current phase and another 
phase. In order to find out whether the latter is the low- 
density phase or the high-density phase, we studied the 
trend of variation of the density profile across the phase 
boundary (see figO. If the average density increases 
and, finally, saturates in the maximal current phase (as 
observed in the inset of figlHta)) while the current reaches 
its maximum value, we identify it with the transition 
from the low-density phase to the maximal current phase 
(as in fig llip . On the other hand, gradual decrease of the 
average density and its eventual saturation (as observed 
in the inset of figlUb)), while the current approaches its 
maximum value, indicates a transition from the high- 
density phase to the maximal current phase (as in fig |12[) . 

Although the extremum current hypothesis [62[ is be- 
lieved to be exact, at least in the limit L — > oo, our 
results on the phase boundaries are approximate because 
we have used the mean-field estimates (fll?|) and (|33p for 
p* and p_, respectively, in equation (|34[) . 



VII. SUMMARY AND CONCLUSION 

TASEP is the simplest model of systems of interact- 
ing "self-propelled" particles. Interestingly, the TASEP 
itself was developed originally [l(| to describe traffic-like 
collective movement of ribosomes on an mRNA strand. 
The physical properties of the TASEP and similar driven- 
diffusive lattice gas models have been investigated exten- 
sively in the recent years using the techniques of non- 
equilibrium statistical mechanics [13, [5i§]- Success of 
TASEP-like models in vehicular traffic Q has not only 
led to the recent modelling of molecular motor traffic by 
suitable extensions of TASEP @,!,!,mmip,[6|L but 
has also revived interest in ribosome traffic jlOl. Illl [H, 

mm mm Hi- 

In reality, a ribosome is not just a "particle" but one 
of the most complex natural nanomotors [l[ . An under- 
lying implicit assumption of the TASEP type models of 
ribosome traffic is that the numerical value of the hop- 
ping rate q is determined by the main rate- limiting step in 
the mechano-chemical cycle of a ribosome. In contrast, in 
this paper we have developed a model of ribosome traffic 
during protein synthesis by explicitly incorporating all 
the major steps in the mechano-chemical cycle of each 
ribosome, in addition to the mutual exclusion of the ri- 
bosomes arising from their steric interactions. Thus, our 
work can be viewed as an interesting biologically moti- 
vated extension of TASEP to an exclusion process for 
extended particles with "internal states" . Exclusion pro- 
cesses with "internal states" have begun receiving atten- 
tion in very recent literature [6!| [7(| • 

We have calculated, both analytically and numerically, 
the flux of the ribosomes, which is directly related to the 
rate of protein synthesis. We have investigated how the 
rate of protein synthesis is affected by the variation of 
the rate constants for the various steps of the mechano- 
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chemical cycle of individual ribosomes. We have demon- 
strated that, with the increase of the numerical value of a 
rate constant, the current eventually saturates when the 
corresponding step of the mechano-chemical cycle is no 
longer rate limiting. 

We have also calculated the average density profiles of 
the ribosomes on the mRNA track in all the dynamical 
phases of the system. Using a few real mRNA sequences 
for E-coli, we have demonstrated that, because of the se- 
quence inhomogeneity, the rate of protein synthesis from 
real mRNA templates is slower than that from a hypo- 
thetical homogeneous template. Besides, we have calcu- 
lated the flux in real time (unlike arbitrary units used in 
most of the earlier works). Our theoretical estimates for 
the rates of protein synthesis are in good agreement with 
the corresponding experimental data. 

Our work also elucidates the nature of boundary- 
induced non-equilibrium phase transitions in a 
biologically-motivated driven-diffusivc lattice gas model 
[27| . We have determined the phase boundaries on the 
phase diagrams for our model by using the extremum 
current hypothesis 62]. Following the traditional ap- 
proach to phase diagrams of TASEP under OBC, all the 
earlier works on TASEP like models of ribosome traffic 
reported phase diagrams in the a — (3 plane. But, we 
have plotted the phase diagrams in planes spanned by 
experimentally accessible parameters that include the 
concentrations of aa-tRNA and GTP-bound elongation 
factors. We hope the predictions of our theoretical 
model will stimulate further experimental studies for 
more accurate quantitative data. 
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Solving (|35|) gives 



[A]o - [A] 
[A]o 



-kAt 



(36) 



where [A] Q is the concentration of A at t = 0. The left 
hand side of (|36|) gives the fraction of A molecules re- 
acted in time At, and is thus the probability that a single 
molecule of A will be transformed into B, in time At. If 
this time interval At is very small, we can expand the 
right hand side of equation (|3"6"1) . Differentiation of this 
with restect to time gives the probability of transition 
per unit time: 



d 

— lim(l 

dt t^o y 



-kt\ 



(37) 



If Pa is the probability of finding the molecule in state 
A, then the final master equation, according to (|3T[) is 



BPa 
dt 



-u) A ^ b Pa 



(38) 



Appendix B: Results in the special case t = 1 

Consider the special case I = 1 of our model of ribo- 
some traffic (with 5 "internal" states for each ribosome) 
on a mRNA of L codons. This model is not equivalent to 
the TASEP-like model of Lakatos et al. 7l| where parti- 
cles (without internal states) hop on a lattice of 5L sites 
with spatially periodic hopping rates. Under periodic 
boundary conditions, the factor Q(i\i + £) in equations (j4|) 
and © reduce to the simpler form Q(i\i+1) = 1— P(i+1). 
Thus, for I = 1, in the steady-state with PBC, 



Pr, 



P 



i + n h2 (i-P) 



(39) 



and. hence, 



Appendix A: Rate constant versus probability 

Consider a chemical reaction: A — >B, with a rate con- 
stant k. Thus, 



d[B] 
dt 



= k[A] = 



d[A] 
dt 



(35) 



J = LO h2 P 5 (l - P) = 



Uh.2p(± ~ P) 

i + n h2 (i- P ) 



(40) 



In the special case k e ff — > 00, where u>h2 = q is non- 
zero and finite, Qh2 —* 0; in that case the expression ([4*0]) 
reduces to the corresponding formula @ for TASEP. 
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